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FOREWORD 


The present report is one of a series of six reports, published simul- 
taneously, which describe analyses and computational procedures for: 1) pre- 

diction of the in-depth response of charring ablation materials, based on one- 
dimensional thermal streamtubes of arbitrary cross-section and considering 
general surface chemical and energy . balances, and 2) nonsimilar solution of 
chemically reacting laminar boundary layers, with an approximate formulation 
for unequal diffusion and thermal diffusion coefficients for all species and 
with a general approach to the thermochemical solution of mixed equilibrium- 
nonequilibrium homogeneous or heterogeneous systems. Part I serves as a 
summary report and describes a procedure for coupling the charring ablator 
and boundary layer routines. The charring ablator procedure is described in 
Part II, whereas the fluid-mechanical aspects of the boundary layer and the 
boundary-layer solution procedure are treated in Part III. The approximation 
for multicomponent transport properties and the thermochemistry model are 
described in Parts IV and V, respectively. Finally, in Part VI an analysis 
is presented for the in-depth response of charring materials taking into ac- 
count char-density buildup near the surface due to coking reactions in depth. 

The titles in the series are: 


Part I 

Part II 

Part III 
Part IV 

Part V 

Part VI 
This 


Summary Report: An Analysis of the Couple'd Chemically Reacting 

Boundary Layer and Charring Ablator, by R. M. Kendall, E. P. 
Bartlett, R. A. Rindal, and C. B. Moyer. 

Finite Difference Solution for the In-depth Response of Charring 
Materials Considering Surface Chemical and Energy Balances, by 
C. B. Moyer and R. A. Rindal. 

Nonsimilar Solution of the Multicomponent Laminar Boundary Layer 
by an Integral Matrix Method, by E. P. Bartlett and R. M. Kendall. 

A Unified Approximation for Mixture Transport properties for Multi- 
component Boundary-Layer Applications, by E. P. Bartlett, R. M. 
Kendall, and R. A. Rindal. 

A General Approach to the Thermochemical Solution of Mixed Equilib- 
rium-Nonequilibrium, Homogeneous or Heterogeneous Systems, by 
R. M. Kendall. 

An Approach for Characterizing Charring Ablator Response with In- 
depth Coking Reactions, by R. A. Rindal. 
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ABSTRACT 


The analysis of a charring ablation material which may undergo subsurface 
coking of the pyrolysis gas is considered. Coking reactions considered include 
thermal cracking of gaseous hydrocarbons resulting in precipitation of carbon 
onto the char layer and, at higher temperatures, the subsequent internal chem- 
ical erosion of the char structure by the gaseous pyrolysis products. A gen- 
eralized type of ablation material is defined which consists of inert, carbon, 
and reactive constituents. Three types of subsurface reactions are considered, 
1) decomposition of up to three organic constituents to form initial char and 
pyrolysis gas products, 2) kinetically controlled decomposition of the pyrol- 
ysis gas resulting in carbon precipitation and char densif ication, and 3) at 
high temperatures, chemical erosion of the subsurface char matrix by the py- 
rolysis products according to the dictates of chemical equilibrium. A model 
is considered for evaluating the pressure distribution through the char layer 
of variable permeability. Differential equations are developed to represent 
the transfer of mass, energy, and momentum within the framework of the postu- 
lated phenomenological model, and the differential equations are subsequently 
cast into finite difference form suitable for coding into a computer program. 
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AN APPROACH FOR CHARACTERIZING CHARRING ABLATOR RESPONSE 
WITH IN-DEPTH COKING REACTIONS 

SECTION 1 
INTRODUCTION 


The theoretical characterization of ablation phenomena has received 
considerable attention during the past decade. Part II of the present series 
(Ref. 1) includes a partial bibliography on the subject which includes 
reference to a number of investigations directed toward mathematical char- 
acterization of various aspects relating to the subsurface behavior of 
particular material classes, to thermochemical interactions between a charring 
ablator and its environment, and to various material removal regimes which 
dictate the magnitude of surface recession. This report is directed toward a 
potentially important aspect of charring material response that has received 
little attention previously. Consideration is given to the ablation of a 
charring material which may undergo in-depth coking reactions and a mathematical 
formulation is proposed for systematically evaluating the extent and effects of 
forward and reverse coking reactions within the constraint of conserving chemical 
elements in-depth and at the heated surface. The equations are cast in finite 
difference form appropriate for coding into a computer program. 

The mathematical model proposed herein represents an extension of that 
presented in References (1-4) in that certain chemical reactions between gaseous 
organic pyrolysis products and the char layer are considered. Experimental 
results presented in References 5 and 6 reveal that significant densif ication 
of the char layer may occur near the heated surface. The model described in 
this report is based upon the assumption that the densif ication results from 
a certain class of chemical reactions which cause a transfer of carbon between 
the pyrolysis products and the char layer. These reactions include the pre- 
cipitation of carbon from the hydrocarbon containing pyrolysis products with 
attendant deposition upon the char (coking) , and the reverse reaction evidenced 
by erosion of the carbonaceous char with attendant addition of carbon to the 
gaseous pyrolysis products. The forward and reverse coking reactions may occur 
in the low and high temperature regions of the char layer respectively, and 
are of interest because the permeability of the char layer is decreased by 
coking which may result in high gas pressures in depth. High pressure in depth 
may give rise to excessive char stress which may produce catastrophic failure 
of the char layer. The technique described herein includes consideration of 
internal pressure build-up resulting from pyrolysis product flow through a 
char layer, of variable permeability. 



The model selected to represent the material state in depth is presented 
first, in Section 2; the conservation equations for chemical species, energy, 
and momentum are presented next, in Section 3; and are followed, in Section 4, 
by the finite difference representation of the differential equations. 

SECTION 2 

DESCRIPTION OF THE PHYSICAL PROCESS 

In this section the type of ablative material being considered is 
described and the general types of reactions which are allowed to occur are 
identified. The model selected is believed to enable consideration of the 
most important thermochemical interactions controlling the in-depth response 
of a wide variety of ablation materials currently being considered for heat 
shielding applications. 

In its undecomposed state the ablation material is taken to be composed 
of two basic types of constituents: 1) inert, and 2) reactive. The inert 
constituents will consist of materials which are not permitted to undergo 
molecular changes in depth, e.g., silica or other metal oxide reinforcements. 

The reactive constituents may consist of organic materials, carbon or graphite 
reinforcements, and water of crystalization of reinforcing fibers, for example. 
Carbon and graphite are included in the list of reactive constituents because 
they may be vaporized in depth or be eroded chemically by the gaseous products 
of other reactive constituent pyrolysis products. The following, idealized, 
irreversible reaction characterizes the initial decomposition of the composite. 

Inert + Reactive -► Inert + Carbon (S) + Gas (0) (1) 

As noted from this relation, the inert constituent does not take part in the 
reaction, but is simply transported from a constituent in the virgin plastic 
to a constituent in the initial decomposed material. The reactive constituents, 
on the other hand, do undergo a change in molecular configuration and phase, 
however, the products of this initial reaction may consist of only two 
constituents, solid carbon, and an initial pyrolysis gas ( gas (0)). Reaction 
(1) should be looked upon as a reaction which splits the virgin material into 
three distinct parts, each having a fixed quantity of chemical elements. 

Previous investigations (e.g. , Refs. 1-4) consider no further chemical 
reactions between the pyrolysis products and the char. Indeed, further 
decomposition of the pyrolysis gas to yield a different gas composition may 
be treated by previous models, but decomposition to yield precipitation of 
carbon onto the char (coking) may not be treated with these models. In the 
present formiilation , the initial pyrolysis is treated by Reaction (1) . The 
gas elemental composition resulting from this initial decomposition is fixed 
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but further reactions are allowed and will be discussed subsequently. , It 
should be pointed out that even though the initial gas composition is ‘fixed 
with respect to its elemental composition, its molecular composition will 
depend upon the temperature and pressure at the point where decomposition 
occurs. Ideally, the elemental composition of the gas evolved at a point 
should depend upon the particular reactive constituent which is undergoing 
decomposition at that point, e.g. for nylon- phenolic ; nylon would decompose 
first to yield something like C^H^ON , which would be followed by decom- 
position of phenolic to yield something like C-jH^O at higher temperature. 

The inclusion of such detailed considerations in a computational scheme, 
however, would represent a rather extensive effort which is not believed 
warranted for the small gain in accuracy which would result. The approach 
taken, therefore, is to consider the elemental composition of the initial 
pyrolysis gases to represent the elemental composition of all gaseous pyrolysis 
products taken collectively. The initial off-gas elemental composition, 

Kk g o ’ ma Y then be obtained by subtracting the quantity of chemical elements 
contained in a laboratory produced char from the chemical elements contained 
in the virgin plastic. 


Vg.c 


= m p K ^P 

m p ' m c 


( 2 ) 


(k = H,C,N,0, e.g.) 


where 


K. 


k , g , o 


m 


P 


m 


c 


kp 


x kc 


represents the mass of element k per unit mass of gas initially 
evolved 

mass of virgin plastic sample 

mass of char after decomposing in an inert atmosphere at a moderate 
temperature (2,000 - 3,000°R) 

mass of element k per unit mass of virgin plastic 
mass of element k per unit mass of char. 


After the initial decomposition gas is formed, it will percolate through 
the char layer toward the heated surface which will result in an increased gas 
temperature and decreased pressure. The change in pressure and temperature 
will cause the initial gas products (gas (0)) to undergo numerous chemical 
reactions as they pass through the char. The reactions considered here fall 
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into three general categories: 


1. Decomposition of the gas including thermal decomposition of 
high-molecular-weight hydrocarbons and dissociation of CC^ > 
f^O, and H 2 for example. 

2. Further decomposition of the hydrocarbons resulting in 
precipitation of carbon (coking) on the adjacent char 
passages resulting in char density buildup. 

3. Chemical erosion of the char layer (below the heated surface) 
by the gases including carbon vaporization resulting in a 
char density reduction near the heated surface. 

The three reaction regimes are represented on the following sketch for decom- 
position of a hypothetical material. 

Char 

Erosion ... . 


Q_ 

>1 

4J 

W 

C 

d) 

Q 


0 

Distance below heated surface 



The classification of reactions above corresponds to the order in which the 
various types of reactions would be expected to occur as the gas passes 
through the char. The first class of reactions may be looked upon simply as 
the gas reacting with itself so no change in the concentration of chemical 
elements in the gas results, K^g 58 ^go * The seconc ^ two classes of reactions, 
however, will result in a transfer of carbon elements between the gas and the 
porous char. With regard to the elemental composition change these reactions 
may be considered reversible. 


gas (0) * — - gas + carbon (S) (3) 

3 


4 



where the forward reaction corresponds to moderate temperature, type 2 reactions, 
and results in a precipitation of carbon from the initial pyrolysis gas. The 
subsequent, type 3, high temperature reactions, result in char erosion with 
attendant addition of carbon to the gas. 

An understanding of the detailed kinetic mechanisms required to char- 
acterize these reactions is not presently in hand; however, some qualitative 
information is available upon which a crude model may be formulated. The 
specific information relating to each reaction type is presented briefly and 
the physical model adapted to characterize each reaction regime is described 
in the following paragraphs. 

Type 1 Reactions - If the subsurface composition is computed on the 
basis of chemical equilibrium considerations a far more dense char is predicted 
to occur than is observed from char density measurements (Ref. 2) . It may be 
concluded either that condensed phase carbon is formed but does not stick to 
the char, or that the high molecular weight hydrocarbons do not decompose 
according to the dictates of chemical equilibrium. The latter possibility 
seems more probable but conclusive experimental evidence on this matter is 
lacking. It would seem appropriate to compute the initial gas phase elemental 
composition from Equation 2, and to evaluate the gas phase molecular composi- 
tion by assuming gas phase chemical equilibrium and not consider condensed 
phase carbon as a possible product. Additional constraints must be placed 
upon the equilibrium calculation in order to eliminate other unlikely species 
from occurring as well. For example, mass spectrometer measurements presented 
in Reference 7 do not indicate the presence of any C^qH, , or CgFi, whereas 

these species are predicted to occur to a significant extent from equilibrium 
calculations for the resin system which was analyzed. By employing experi- 
mental techniques such as those reported in Reference 7, one should be able 
to identify those molecular species which will occur. The distribution of 
these experimentally verified species may then be obtained from equilibrium 
considerations . 

Type 2 Reactions - A certain amount of experimental evidence exists 
(Refs. 5 and 6) which indicates that char densif ication may occur between the 
organic decomposition zone and the heated surface. It appears that this char 
densif ication is a result of deposition of condensed phase carbon from the 
hydrocarbons in the organic pyrolysis gas products. This seems reasonable 
since, as indicated above, the gases contain far more carbon than would exist 
if equilibrium were achieved. As the gases approach the heated surface their 
temperature is increased and the rate at which equilibrium is approached 
increases. In the present study, the rate at which equilibrium is approached 
is represented by a kinetic equation of the Arrhenius form. The rate at which 
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coking reactions may proceed is taken to be proportional to a forward rate 
coefficient and a driving potential represented by the difference in chemical 
composition between the actual gas composition and the equilibrium gas com- 
position. Consider the reaction rate of coking methane gas, for example 1 


coke 


CH„ 




(4) 


where K 

P 


is the equilibrium constant for the coking reaction 


ch 4 H 2h 2 + C (S) 


Expressing the forward rate coefficient in Arrhenius form yields: 


k 




-E./R T 

A . e 1 
1 


(5) 


( 6 ) 


where the frequency factor (A.) and activation energy (E.) are characteristic 

th ^ ^ 

of the l reaction such as Reaction 5. The total coking rate from the gaseous 

pyrolysis products would be given by the Siam of the coking rates attributable 

to each reaction. 

m coke ” 

i 

where i includes coking reactions for all hydrocarbon species present. The 
above technique for evaluating the coking rate, aside from requiring time- 
consuming computation procedures, would require specification of the physical 
constants, and E^ , for each of the many reactions being considered. The 

rate would also depend critically upon the specific concentrations of each 
molecular species, and, as indicated earlier, the methods for their evaluation 
are not presently in hand. Alternatively, an approximate procedure is pro- 
posed here which has the following features: 

1) It is simple enough to be included practically in a 
charring ablation solution. 


I 


m 


coke . 


(7) 


The reaction is presumed ideal so the partial pressure exponents are equal 
to the stoichiometric coefficients. 
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2} It approaches the coking rates which would be predicted 
by the more complete equations above when kinetics are 
relatively slow or fast. 

3) It is based upon parameters which may be controlled in 
a laboratory experiment to derive data on the coking 
process . 


The present procedure consists of specifying the coking rate as the product 
of a single forward rate coefficient and a carbon mass fraction "coking po- 
tential" which represents the net effect of all coking reactions such as 
Equation (6) . 


m 


coke 


k f^ K cg K cgE^ 


( 8 ) 


where the forward rate coefficient is expressed in Arrhenius form as above, 
and the driving potential is represented by the difference between the 
elemental carbon mass fraction of the gas and that which would exist if 
equilibrium were achieved. 


Figure 1 depicts the elemental carbon content of a typical pyrolysis gas 
for the following situations: 1) for the initial pyrolysis gas (Equation (2) , 

2) for the pyrolysis gas if equilibrium were achieved, and 3) for a typical 
case utilizing the above model, all as functions of temperature. It is noted 
that at low temperatures (regime 1) the elemental carbon mass fraction remains 
unchanged until regime 2 is reached, at which point kinetically controlled 
coking reactions will cause the elemental carbon content to decrease until the 
equilibrium composition is reached. This marks the onset of regime 3 which is 
characterized by addition of carbon to the gas from the char. 

Type 3 Reactions - In the event local char layer temperatures much in 
excess of 4500°R are achieved, it is probable that chemical equilibrium will 
be achieved between the pyrolysis gas and the char, in which case the coking 
potential will reach zero. For this reason type 3 reactions are presumed to 
occur in chemical equilibrium and sufficient carbon will be added to the gas 
from the char to maintain this equilibrium. This char erosion will result in 
a char density decrease which may, in an extreme case, cause the char density 
to approach zero. 

SECTION 3 

FORMULATION OF THE SUBSURFACE CONSERVATION EQUATIONS 

The differential equations which characterize the conservation of species, 
momentum , and energy within a charring ablation material are presented in this 
section. The equations are based upon the physical model described above. 
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3.1 SPECIES CONSERVATION 

In order to properly characterize the rates of plastic decomposition 
and pyrolysis gas coking in depth it is necessary to accurately assess the 
state of both the ablation material and the pyrolysis gases. The state of 
the ablation material depends upon the relative quantities of each constituent 
(reactive, inert, and carbon) and temperature whereas the state of the 
pyrolysis gas depends upon the gas carbon content, temperature, and pressure. 
Temperatures are obtained from the energy equation and pressures from the 
momentum equation, both of which are presented subsequently. The species 
conservation equations, presented here, relate the rate of change of ablation 
material composition and pyrolysis gas composition to decomposition and coking 
events throughout the decomposition zone and char layer. 

As indicated earlier (Reactions 1 and 3) the ablation material may, at 
any point in space and time, be represented by a mixture of inert, reactive, 
and carbon species. The geometrical configuration considered is shown in 
Figure 2. Considering a control volume of extent Ady, the total mass may be 
expressed as 


m 


m T + m r + m c + m g 


(9) 


where 


mass inert constituents 
mass of condensed carbon 


mass of gas 


z 


m 


i ’ 


mass of all 


i reactive constituents 


It is noted that the mass of gas, m g in the control volume is, in general, 
much smaller than the other masses so the total mass in the control volume 
may be written 


m “ 


E 


m 


a 


( 10 ) 


where a - I, c, r 
and r - 1, 2, 3 

It will be convenient later to consider a coordinate system fixed to the 
receding surface (x m const) . With this in mind, the mass conservation 
equation is written for the moving coordinate system by differentiating equa- 
tion (10) at constant x . 
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■ M ? ■*)* 


(ID 


The mass change rate of a constituent, a , may result from chemical reactions 
(organic decomposition, coking, or char erosion) and convection resulting from 
coordinate system movement. The mass change rate at constant x may be 
expressed in terms of the change resulting from chemical reactions and the 
change resulting from coordinate system motion. Functionally, the mass of 
constituent a may be written in terms of time and position. 


m a * m a (0 ’ y) 


Applying the chain rule yields 


dm a 



d 0 + 



dy 


Differentiating with respect to time at constant x obtains 



( 12 ) 


The second term may be related to the surface recession rate by considering 
the coordinate system dependence of x on y . 


y * x + S 


(13) 


From which, 



but the surface recession, S , may be considered as a function of time alone, 
so: 


bs 

Je 



(14) 


Substitution of Equation (14) into Equation (12) and noting that differentia- 
tion with respect to x or y at constant time is equivalent, yields a 
conservation equation for condensed phase species below the surface. 
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s 


(15) 


Je~) 


X 


dm \ din 

^r) y + ssr 


e 


The conservation equation for gaseous species in depth follows directly 
from consideration of the gas generation rate resulting from decomposition and 
the gas-to-condensed-phase transfer associated with coking reactions. 



e 


d 

Te 


<Pa A >y 


(16) 


As indicated in Section 1 (Reaction 1) the reactive constituents are presumed 
to decompose into a gas having a fixed elemental composition, given by Equa- 
tion (2) , and that subsequent gas phase reactions (coking or erosion) will 
have the sole effect of removing or adding elemental carbon to the gas. It 
is therefore necessary to keep track of only one element in the gas phase, 
elemental carbon. The gas phase conservation equation for elemental carbon, 
a - c , is 


^y ( *cgV e " fs ( K C P A) y 


(17) 


The overall mass conservation equations for condensed and gas phase species 
follow directly from summing Equations (15) and (16) over all "a” species. 


dm \ m dm \ dm 

H ®|y H 


(18) 


and 


dm 

■Sy 1 


<P A >y 


(19) 


dm. 


Evaluation of the convection term, s > in the species conservation equa- 

tion is relatively straightforward; however, the term representing species 

Sm a \ 

change rates associated with chemical reactions, -gg— , requires rather 

extensive consideration. For the inert species, a ■ I , the mass change rate 
resulting from chemical reactions is zero. 


bm i\ 

■Se~j 


y 


o 


(20) 
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The treatment given the other constituents, a * r and c , is described in 
the following two subsections. 


3.1.1 Species Conservation for Reactive Constituents (a » r) 

The reactive constituents, r , are defined here as those constituents 
which may decompose to form gaseous and/or carbonaceous pyrolysis products 

Reactive -► Carbon (S) + Gas (0) (21) 


The above, irreversible reaction is taken to proceed at a rate governed by an 
equation of the Arrhenius form for each of up to three reactive constituents 


dm r \ V ' dK? \ 

I • t 9 m io d0 1 

I y 1 1 > 2 , 3 / y 


where 




Ho- 


ly 


-k. 


-E^/R T 


(kJ) n± 


(23) 


and K? 



mass of constituent i per unit mass of constituent 
i prior to decomposition. 


The form selected to represent the decomposition rate arises from a desire to 
express the kinetic coefficients in terms of quantities readily measureable 
in the laboratory. Equation (23) is developed in Appendix A in terms of 
quantities readily derived from TGA (thermogravimetr ic analysis) data. 

Equations (22) and (23) enable evaluation of the term in the 

species conservation Equation (15) for decomposable reactive constituents 
(a » r) . in order to evaluate this term for carbon (a * c) it is necessary 
to first consider all reactions which may involve carbon. This includes both 
the decomposition reactions just described, and coking reactions (forward and 
reverse) . Coking reactions are discussed next. 


3.1.2 S pecies Conservation for Carbon (a a c) 

Carbon is generated in the charring material matrix when the virgin 
material decomposes as described above, and when coking reactions cause a 
transfer of carbon from the pyrolysis gas to the char (zone 2 described in 
Section 2 above) . Carbon is removed from the charring material when reverse 


11 



coking reactions effect a transfer of carbon from the char to the pyrolysis 
gas (zone 3) . The conservation of carbon species associated with decomposition 
is considered first, in Section 3. 1.2.1, and is followed, in Section 3. 1.2.2, 
by carbon conservation relations associated with forward and reverse coking 
reactions . 


3. 1.2.1 Organic Decomposition 

From Equations (A-3) and (A-6) developed in Appendix A, it may be shown 
that the rate of production of carbonaceous residue from the i^ decomposition 
reaction is related to the rate of decomposition of constituent i , in the 
following manner 


dm 


"Se" 


Cl 



*Hy 


(24) 


where K° = — is the mass of carbonaceous residue resulting from com- 

C if ra io 

plete decomposition of constituent i , per unit mass of constituent i prior 
to decomposition. Summation of Equation (24) over all i yields the carbon 
production rate as a result of organic decomposition reactions. 



(25) 


where the subscript d signifies carbon production associated with organic 
decomposition (Reaction 21) . Carbon addition to and removal from the char 
(coking reactions) are considered next. 


3. 1.2. 2 Coking Reactions 

Coking reactions, as employed here, refer to all chemical reactions 
which result in a transfer of carbon atoms between the pyrolysis gases and 
the char structure. This includes both type 2 and type 3 reactions introduced 
earlier in Section 2. Type 2 reactions include the normal class of coking 
reactions characterized by the deposition of carbon upon the char resulting 
from hydrocarbon decomposition. Type 3 reactions occur at much higher tem- 
perature where equilibrium dictates a reverse reaction characterized by 
erosion of the char resulting from vaporization and chemical attack by the 
pyrolysis gases. 
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3. 1.2. 2.1 Forward coking reactions 

The physical model to be employed for characterizing the rate at which 
carbon is coked from the gas was introduced above, in Section 2. It is 
assumed that the pyrolysis gas composition changes at a rate dictated by an 
equation of the Arrhenius form. 



-E c /R 


(P) C (K 


eg 


K cgE^ 


(26) 


where 



K cg E 


mass fraction of carbon in the pyrolysis gas 
independent of molecular configuration 
pressure 

kinetic coefficients for coking reactions 

mass fraction of carbon in the pyrolysis gas 
which would exist if equilibrium were achieved 
at the local temperature and pressure 


Equation (26) expresses the rate of change of carbon content in the gas in 
terms of a driving potential which approaches zero as equilibrium is achieved. 
The driving potential is noted to depend in a direct manner upon the system 
pressure as is believed appropriate for gas phase reactions. Partial 
differentiation of the gas carbon content is performed at constant y and 
constant (composition of the plastic-char matrix where coking is 

occurring) . It is expressed in this manner because coking is allowed to occur 
in regions of the char still undergoing dec apposition . Since local pyrolysis 
results in the addition of fresh gas (having a carbon content, K cgQ ) to the 
pyrolysis gas stream, it too affects the gas composition change rate. 

Equation (26) reflects the gas composition change rate resulting from coking 
alone. The total carbon change rate of the pyrolysis gas resulting from both 
decomposition and coking is given by adding the effects of that due to 
decomposition (Equation (22) and coking (Equation (26)). This addition is 
given consideration subsequently, in the finite difference formulation 
(Section 4) . 


3. 1.2. 2.2 Reverse coking reactions 

Referring to Equation (26) , it is noted that the coking rate approaches 
zero as gas phase equilibrium is achieved, i.e., as ^ K cg "* K C gE^ 0 * When 
the coking potential reaches zero, (K cg - K cgE^ = 0 » 9 as :i ' s assumed to 
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be in equilibrium with the char layer: 


K cg K cgE 


(27) 


where the pyrolysis gas equilibrium carbon content may be expressed as a 
function of temperature and pressure alone. The rate of change of elemental 
carbon content in the pyrolysis gas is then related to the temperature and 
pressure gradient. 




(28) 


The rate of char erosion by the pyrolysis gas (reverse coking) may then be 
obtained from Equations (28) and (17). This matter is given further con- 
sideration subsequently, in Section 4. 


3.2 MOMENTUM EQUATION 

The momentum equation is considered in depth in order to relate the gas 
flow rate in the char layer to the pressure distribution throughout the char 
layer. It is desired to evaluate the pressure at all points in the char 
layer for two reasons: 1) the coking reactions, both forward and reverse, are 

pressure dependent, and 2) evaluation of the pressure distribution through 
the char layer will enable determination of char stresses in an approximate 
manner . 

In order to obtain an expression relating local pressure in the char 
layer to other pertinent variables it is useful to examine experimental data 
taken for a range of variation of the pertinent parameters of interest, and 
to deduce the most significant effects by generalizing this data. 

The data presented by Green (Ref. 8) is employed for this purpose. 

Green presents a compilation of data for flow of various gases through a wide 
variety of porous media for a large range of flow conditions. In order to 
relate pressure gradients to viscous and inertial forces he employs the 
correlation equation first proposed by Reynolds. 


dp 

dy 


2 

a (rv + ppv 


(29) 


where the first term represents viscous forces and the second, momentum forces. 
The quantities a and 3 are empirical coefficients, is the gas viscosity, 
and p is the gas density. The velocity in Equation (29) is a "superficial 
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velocity ' 1 defined on the basis of the gas flow rate per unit projected area 
in a plane normal to the velocity vector. 


v = 



(30) 


Referring to Equation (29) the ratio of inertial to viscous forces may be 
written as: 


Inertial Force _ p g 

Viscous Force apA ' 

where m g = pv is the gas flow rate per unit total area. The compilation of 
data presented by Green includes tabulations of the empirical coefficients 
a and 0 for a wide variety of porous media including packed beds of 
irregular and spherical particles ranging in nominal size from 0.08 inch to 
0.1875 inch and for close packed wire screens ranging from 60 to 5 mesh. The 
nominal range of porosities included varies from 0.3 to 0.8. The experimentally 
derived coefficients a and 0 are shown in Figure 3 along with a straight 
line fit to the data. 


a = 0.794 x 10 P 
—2 —1 

where a has units of ft and 0 has units of ft . The correlation is 
not excellent, but appears appropriate for order-of-magnitude considerations. 
Substituting this relation into Equation (31) obtains 


Inertial Force 
Viscous Force 


1.26 x IQ - * 


UA 


For the temperature range of interest (2000°R - 5000°R) the gas viscosity 

_4 _4 

will range from about 0.3 x 10 to 0.5 x 10 lb/ft-sec. Considering a gas 
viscosity of 0.4 x 10 lb/ft-sec, the ratio of inertial to viscous forces 
may be written: 


Inertial Force 3 ^2. 

Viscous Force ' A 

• 2 

where m^/A has the units of lb/ft -sec. It may therefore be concluded that 
for pyrolysis gas flow rates of 0.01 lb/ft -sec and less, the inertial terms 
may be ignored in which case Equation (29) reduces to Darcy's Law. 
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dp 

= 


where a -1 is the permeability. It is recommended, however, that the more 
complete correlation equation (29) be employed for pressure drop calculations 
in the char layer of ablating materials since it is valid for a wider range 
of conditions of practical interest. 


3.3 ENERGY EQUATION 

The subsurface energy equation for a charring ablator may be written in 
terms of a stationary coordinate system (see Fig. 2) . 


storage 


conduction 


convection 




a 

^7 


XA 


aT 

^7 


a 

"57 


( "gVe 


(32) 


In Equation (32) the rate of energy storage term is expressed in terms of the 
material density, p , and specific enthalpy, h , at a point. The implicit 
assumption in such a formulation is that the density and enthalpy are sufficient 
properties to represent the material state. This is indeed possible when 
considering simple irreversible decomposition of the plastic into char and 
gas (see e.g. Refs. 1-4) , but when further reactions accompanied by mass trans- 
fer between the char and gas are considered, the density and enthalpy alone 
are not sufficient properties to establish the material state. It is there- 
fore desired for the present analysis to express Equation (32) in terms of 
quantities sufficient to define the state for the situation where mass trans- 
fer may occur between the pyrolysis gas and the char. It is also most 
expedient if these quantities are readily accessible in the computation 
scheme. These considerations have led to the selection of the mass of 

constituent a , in a control volume (m ) and specific enthalpy of component 

a 

a (h ) as the most convenient state properties. 

3 

Equation (32) may be rewritten in terms of these properties by integrating 
over the extent of the control volume shown in Figure 2 


E 


Te 


(m h ) 

N a a' y 


XA ' 
kA 

V ^1 


^“gVe 


(33) 


where a = I, c, or r , and the differences on the right hand side are between 
the top and bottom of the control volume. In order to develop a practical 
numerical solution it is convenient to consider a coordinate system fixed to 
the receding surface (a moving control volume) . For this purpose it is 
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desired to transform Equation (33) , which is written for a control volume, 
y = constant , to an equation written for a moving coordinate system, 
x - constant. The storage term on the left side of Equation (33) may be 
related to its counterpart in the moving coordinate system in the following 
manner. It is noted that the mass of constituent a and its enthalpy may 
be expressed purely as functions of space and time. 

m a h a = f(y,e) 


so: 


"S^ (m a h a ) x = l^ (m a h a ) y + 0 l$) x (34) 

The x and y coordinates and their derivatives are related through Equa- 
tions (13) and (14) 



(35) 


Utilizing this relation and noting that differentiation with respect to x 
or y at constant time (0) is equivalent , enables expressing Equation (34) in 
the following manner: 


a 


(m h , 

a a' y 


x _ a 


(m h ) 
v a a' x 


® b {m 


a h a>0 


Substituting the above into Equation (33) yields a subsurface energy equation 
in terms of the moving coordinate system. 


z 




II 


kA 


dT 


III 


+ ® ax (m a h a ) e 


IV 


(36) 


Terms I and IV involve time and space derivatives of the control volume 
mass and enthalpy. It will be convenient for subsequent finite difference 
formulation to express enthalpy derivatives in terms of temperature deriva- 
tives. Also, the mass derivatives both include cross sectional area deriva- 
tives which may be eliminated. The following manipulations accomplish the 
desired result. 


Term I will be considered first. 


d 

c>e 


(m a h a } x 


m„ 


Sh a\ 

cnr 



(37) 
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The enthalpy for any constituent, a - I, c, or r , may be expressed in terms 
of its chemical and sensible parts 




h, = AH^ + / C dT 


/ 


(38) 


For a given constituent, AH^ = constant, so: 

a 


ah al 


, dT 


(39) 


Utilizing Eouation (35) the mass derivative may be expressed as follows 


dm l dm a \ 

^r) x ^~l y + 


dm ’ 

a 

7e~ 


(40) 


Substituting Equations (39) and (40) into (37) yields the following representa- 
tion for term I in the energy Equation (36) . 


d 


(m a h a>x 


m e 4? 
a p 69 


+ h. 


dm 

rfly ^ S 


•51 


(41) 


Utilizing Equation (39) , term IV in the energy Equation (36) may be 
written as: 


S ^( m A) = S m^ C. 


dT 


a a 


0 


dm\ 


a U P, + s h a 


0 


(42) 


Substitution of Equations (41) and (42) into Equation (36) yields the desired 
form of the subsurface energy equation for a control volume of finite extent 
in terms of a coordinate system fixed to the receding surface. 




a a P a 


S) 8 * 1( W, * 5 5 ’• = P a H) * 5 h * ^)y 


(43) 


It is noted that differentiation of the last term is at constant y. It is 
expressed in this manner for convenience since the decomposition events (a - r) 
are represented at constant y, (see Eq. (22)). 

In the following section a finite difference representation for the 
differential equations presented in this section is proposed. 
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SECTION 4 

FINITE DIFFERENCE FORMULATION 

The need for a computational scheme to mathematically model the thermo- 
chemical response of a charring material experiencing subsurface coking reac- 
tions was introduced in Section 1. In Section 2, a phenomenological model 
was introduced for representing the state of the charring material and its 
pyrolysis products, and for representing the rate at which coking reactions 
may proceed. Section 3 presented abbreviated derivations of the differential 
equations appropriate to representing the transfer of mass, momentum, and 
energy for the chosen phenomenological model. In this section, the equations 
presented above are cast into finite difference form appropriate for coding 
into a computational scheme. 

Examination of the equations presented in Section 3 reveals that finite 
difference formulation of the energy conservation Equation (43) and the pres- 
sure drop relation (Eq.(29)) may be accomplished in a rather straightforward 
manner provided the coefficients appearing in these relations are known . The 
coefficients most usually contain information relating to the state of the 
charring material or the pyrolysis gas products. As the problem has been 
formulated, the state of the charring material is uniquely defined by its 
temperature and composition (i.e., relative quantity of each constituent; 

I, C, or r) , and the pyrolysis gas state is defined by its pressure, tempera- 
ture, and composition (i.e., K^) . The composition of the charring material 
and the gaseous pyrolysis products is obtained from the species conservation 
relations. A finite difference formulation of the species conservation rela- 
tions is presented first, in Section 4.1. Finite difference equations for 
evaluating the char layer pressure drop are presented next, in Section 4.2, 
and are followed, in Section 4.3, by a finite difference formulation of the 
subsurface energy balance. The finite difference representation of the geo- 
metrical configuration results from dividing the material into a number of 
discrete zones or nodes. The generalized coordinate system and nodal iden- 
tification scheme is shown in Figure 4. 

4.1 SPECIES CONSERVATION 

In this section, the differential equations for the rate of change of 
composition of the charring material and the gaseous pyrolysis products are 
applied to a finite control volume in order to obtain algebraic equations 
relating the composition of the material and gas to other problem parameters. 
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The material composition in a control volume is defined by the masses of re- 
active (r) , inert (I) , and carbon (C) in the control volume. The elemental 
composition of the gaseous pyrolysis products is defined by the carbon mass 
fraction (K ) - The material composition in the control volume will be con- 
sidered first, in Section 4.J.1, and will be followed by finite difference 
relations for evaluating the pyrolysis gas flow rate and composition, in Sec- 
tion 4.1.2. 

4.1.1 Material Composition 

The control volume to be considered is shown in Figure 2 and is chosen 
such that its thickness, dx, is constant and its distance below the receding 
(heated) surface is constant (i.e., x = constant) . The composition of the 
material in the control volume may change as a result of either of two events, 
coordinate system movement, or chemical reactions. The total rate of change 
of mass of a constituent in the control volume resulting from these two events 
is given by Equation (15) : 



din \ 

a 

de” 

/ 


9m 


dx 


(15) 


where a = I, C, or r. The first term on the right-hand side represents the 
change in mass of constituent "a" resulting from chemical reactions such as 
decomposition or coking. The second term represents mass change associated 
with coordinate system movement. Finite difference interpretation of the 
chemical reaction term requires that consideration be given to each particular 
class of chemical reactions, i.e., decomposition, coking, or reverse coking. 
Finite difference evaluation of this term will be treated subsequently, in 
Section 4.1. 1.2. The second term, representing coordinate system motion, is 
considered first, in Section 4. 1.1.1. Finally, in Section 4. 1.1. 3, the finite 
difference relations for material composition are summarized. 


4. 1.1.1 Material Composition Change Resulting from Coordinate System 
Motion 

Finite difference evaluation of the convection term (dm^/dx) 0 S is ob- 
tained by considering a control volume of thickness 6^ and cross-sectional 
area A as shown in the sketch below. 
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Surface location at time 0 


0 


Position 
at time 



e 


o + as 


Considering each node to be of homogeneous composition, the mass change of 
constituent "a" in the node "n" during the interval, t\Q , may be obtained by 
evaluating the mass of constituent "a" lost and gained as a result of coordi- 
nate system motion during the time interval. 


mass loss 



AS 

n 


mass gain 



AS 

n+1 


The mass gain rate is obtained by taking the difference of the above and divid- 
ing by the time step, A G. The following representation of the convection term 
results . 


dm 

cl 

dx 



s 


(44) 


It is noted that the subscript “a" refers to reactive (a = r) , inert (a = I) , 
and carbon (a = C) species. in the following section the mass change rate 
associated with chemical reactions is considered. 
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4. 1.1. 2 Material Composition Change Resulting from Chemical Reactions 

Since the inert constituents are not allowed to react, their mass change 
rate due to reactions is zero, as shown by Equation (20) . Finite difference 
representation of the mass change rate resulting from chemical reactions for 
reactive and carbon species may not be zero. They are treated next. 


4. 1.1. 2.1 Reactive constituents (a = r) 

The mass change rate of constituent r resulting from chemical reac- 
tions is given by Equations (22) and (23) . A finite difference evaluation 
for the mass change rate of reactive consitutent "i" may be obtained by inte- 
grating Equation (23) over the time interval A^ 
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-E i /RT 


ae 


where O' indicates evaluation at the time 0 + A0. 

Evaluation of the integral obtains the following for n^ = 1 


K° = K? exp 

In . In 


-E./RT 

- k^e A 0 


(45) 


and, for n^ ^ 1 
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-E./RT 
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i l 
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i -n • 
i 


(46) 


It is noted that the above integrations were performed holding the temperature 
constant over the time interval. Utilization of the temperature at the “old" 
time, 0, results in an evaluation which is implicit in density and explicit 
in temperature. This method has been found to be satisfactory for a rela- 
tively wide range of materials and temperature rise rates (Ref. 1) . Evalua- 
tion of the Equations (45) and (46) may be made fully implicit by utilizing 
the temperature at the end of the time interval, but this introduces rather 
severe computational complications since coupling of the decomposition and 
energy conservation relations is then required. 
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Finite difference evaluation of Equation (23) is given by: 




A0 


Utilizing this relation and Equation (22) enables forming an expression for 
the mass change rate of a reactive constituent, i, in the node "n" . 





K° 

1 


e 


) 


(47) 


where ( is evaluated from Equations (45) and (46) for the reaction order 

n^ = 1.0 and n^ ^ 1.0, respectively. 

Finite difference formulation of the mass change rate of carbon content 
for each node is presented in the next section. 


4.1.1. 2.2 Carbon (a = c) 

The carbon content in a node may change as a result of char formation 
reactions described above, as a result of precipitation of carbon from the 
gas (coking) , or as a result of erosion of the char by the pyrolysis gas (re- 
verse coking) . Finite difference relations for evaluating the change in car- 
bon content of a node resulting from each of these reactions are presented 
in the following three subsections. 


4.1.1. 2. 2.1 Decomposition of reactive constituents 

The rate of change of carbon mass in a node as a result of decomposition 
of reactive constituents is given by Equation (25) 



Decom. 
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(25) 


where subscript Decom. refers to the carbon change as a result of decomposi 
tion alone. Employing the same finite difference representation for the de- 
composition rate as before (Eq. (47)), the carbon mass change rate in the 
node resulting from decomposition is given by: 


^c\ 

in, Decom. 



i=l, 2, 3 




(48) 
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where, as before, K° , is given by Equations (45) and (46) for n. =1.0 

u l 

and ^ 1.0, respectively. The carbon change rate resulting from forward 

coking reactions is considered next. 

4. 1.1. 2. 2. 2 Forward coking reactions 

Equation (26) relates the pyrolysis gas carbon content change rate as a 
result of forward coking reactions to local pressure, temperature, and gas 
composition . 
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(26) 


This equation will be integrated over the extent of the node n, in order 
to relate the gas composition change rate to the carbon precipitation rate 
on the node . 
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where integration is performed over the thickness of the node and t repre- 
sents the transit time for the gas to pass from the bottom to the top of the 
node. The pressure and temperature are treated as constant throughout the 
node so the equilibrium gas composition, K ma Y t>e treated as constant 

over the integration interval as well. For the reaction order different from 
unity (n ^ 1) , integration of Equation (49) yields; 
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The carbon content of the 
by the gas composition of 


gas entering the node, K C gQ> will be represented 
the gas leaving the previous node, n + 1. 
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(51) 
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Similarly, the carbon content of the gas leaving the node, K cg . i s 

identified with the node, n. 


K cg, 0+t 


K 


cgn 


( 52 ) 


The transit time for the gas to pass through the node, t, is related to 
the gas velocity and the nodal thickness 


T 



(53) 


where u g is the actual gas velocity passing through the node and is not to 
be confused with the superficial velocity (v) employed for pressure drop cal- 
culations (Eq. (29)). The actual gas velocity is related to the gas flow 
rate through the node as follows: 


u 


g 


P g eA 


(54) 


where p g is the pyrolysis gas density, € is the interconnected void volume 
porosity, and A is the local cross-sectional area parallel to the heated 
surface. Utilizing Equations (53), (54), and the gas state relation. 
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(55) 


where R = universal gas constant and M = local pyrolysis gas molecular weight. 

g 

Substitution of Equation (55) into (50) results in an expression for the pyrol- 
ysis gas carbon content leaving the node in the presence of forward coking 
reactions . 
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where the subscript coke signifies the gas composition leaving node n as a 
result of coking reactions alone. It is noted that this term must be modi- 
fied if the node is decomposing at the same time it is receiving carbon from 
the gas. Evaluation of the gas composition including all effects is treated 
subsequently, in Section 4.1.2. Equation (56) is valid for a non-unity re- 
action order, n c ^ 1.0. For n c = 1.0, integration of Equation (49) and 
utilization of Equation (55) yields 
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Equations (56) and (57) represent the carbon content of the gas leaving the 
node which would result from forward coking reactions alone for reaction 
orders (n c ) of non-unity and unity, respectively. 

The purpose of this section is to obtain finite difference expressions 
for the carbon content change rate of a node resulting from forward coking 
reactions. This is accomplished by integrating Equation (17) over the node 
and employing Equations (56) and (57) to represent the gas composition. 
Equation (17) relates the gas composition change rate to the char material 
composition change rate. 
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Rewriting Equation (17) yields 
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The first term represents gas composition change as results, for example, from 
coking reactions. The second term represents pyrolysis gas flow rate change 
as results from decomposition of reactive constituents. The latter effect was 
treated above (Eq. (25)). At this point it is most convenient to consider 
the char carbon content change associated with caking reactions alone, and, 
for this purpose the second term may be ignored. Integration of Equation (58) 
then obtains 
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where the integration is performed over the interval from the bottom to the 
top of node n. The integral on the right represents the time rate of change 
of carbon content in the node n as a result of coking reactions. Per f arm- 
ing the integration yields: 
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( 60 ) 


where the subscript ’'coke” signifies the change resulting from coking alone, 

and m is the pyrolysis gas flow rate entering the node n (Jeaving 

9h+l ~ 

node n+1) . Similarly, K cg n+ ]_ represents the composition of t.he gas enter- 
ing node n. The mass fraction of carbon in the gas leaving node n due to 
coking alone, K C g ) co ^ e is evaluated by either Equations (56) or (57) de- 
pending on the reaction order for coking (n^) . 

The mass change rate of carbon in a node resulting from erosion by the 
pyrolysis gas (reverse coking) is considered next. 


4. 1.1. 2. 2. 3 Reverse coking reactions 

The pyrolysis gas composition change between two adjacent nodes result- 
ing from reverse coking reactions is obtained by direct integration of Equa- 
tion (28) which yields: 
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where the equilibrium pyrolysis gas carbon content, K C g E , is evaluated at the 
local temperature and pressure of the designated nodes. Substitution of the 
above into Equation (60) yields the desired expression for reverse coking 
reactions . 
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It is noted that Equation (62) is identical to Equation (60) with the excep- 
tion that the quantity K \ is now replaced with K„_ . it will be 

cg n Jcoke c 9 E n 

most convenient to employ Equation (60) in both events with the interpreta- 
tion that the quantity ^cg n ] co k e represents the elemental carbon mass frac- 
tion in the gas leaving the node n in the presence of coking reactions 
(forward, or reverse) , but not considering the effect of decomposition of 

reactive constituents in the node n. For forward coking reactions, (K 
~ ~ i cg n+l 

K C gE n ) > 0, Equations (56) or (57) are employed to evaluate K cg n ] C oke de P end “ 

ing upon the coking reaction order, and for reverse coking, (^cg n +i ” ^cgE 0/ 
* c 9n)coke = * c 9 E n- 

4. 1.1. 3 Summary of Relations for Evaluating Material Composition 

The finite difference relations developed in this section are summarized 
and combined here. This is accomplished by considering the net change in 
composition of a node resulting from all events taken collectively. Algebraic 
equations for the mass of each constituent at the end of a time step, AO, are 
given in terms of conditions at the beginning of the time step. Inert, reac- 
tive, and carbon constituents are considered first, in Sections 4. 1.1. 3.1 
through 4. 1.1. 3. 3, and are followed in Section 4. 1.1. 3. 4, by relations for 
the total nodal mass. 


4.1. 1.3.1 Inert constituents (a = I) 

Combining Equations (15) and (20) for inert constituents (a = I) yields 



(63) 


Utilizing Equation (44) to represent the right-hand side and expressing in 
finite difference form obtains: 




n+1 


\ 6 n+l 



(64) 


where it is recalled that 6 ' = 0 + . Rearranging yields the desired rela- 

tion for the mass of inert constituent in a node at the end of the time in- 
terval. 


= m T 


n0 ' 


n 0 


n+1 
\ 6 n+1 
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(65) 
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4.1. 1.3.2 Reactive constituents (a = r) 

The mass of reactive (decomposable) constituents in a node at the end 
of a time step depends upon the amount of decomposition and convection result- 
ing from coordinate system motion. Substitution of Equations (44) , (45) and 

(47) into Equation (15) yields the finite difference expression for a reactive 
constituent with reaction order of unity 


n0 1 


m . 
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n 0 


A0 


/ m. 


1 n+l 

K + i 



m . 
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2 K o 

A0 K i 
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1 + exp 


1 -E./RT 
k ± e 1 A0 


( 66 ) 


The mass of constituent i in the node at the end of the time interval is 
then given by: 
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1 + exp k^e 1 A0 


(67) 


The above equation is for reactive constituents having a unity reaction order 
(n^ = 1.0) . Similarly, for constituents having non-unity reaction order. 
Equations (15) , (44) , (46) , and (47) yield: 
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- (1 


-E . /RT 

n^)k_^e 1 A0 


l-Hi 


( 68 ) 


The total mass of reactive constituents in a node is given by: 



i=l ,2,3 


(69) 


4. 1.1. 3. 3 Carbon constituents (a = c) 

The mass of carbon in a node at the end of the time step depends upon 
the magnitudes of convection, decomposition, and coking. Utilizing Equation 
(15) , (44) , (48) , and (60) obtains 
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( 70 ) 



The terms , and ^ C g n ) co ^ e are subject to several interpretations 

follows. For a decomposable constituent, i, having a reaction order of 

K° , is given by Equation (45) 

1 0 • 



K° exp 

x e 


-E./RT 
k^e AO 


and, for non-unity reaction order. Equation (46) is employed. 
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Equations (56) and (57) represent forward coking reactions. Reverse coking 
reactions will occur when (K C g n+ ^ “ K C g En ) _ 0, and for this case 

Kc 9n)coke = K cgE n * 

It is noted that the equations presented thus far in this summary section 
are sufficient for representing the material composition for solution of the 
pressure-drop and energy conservation relations to be presented subsequently. 
Several quantities relative to material composition and density, in addition 
to those presented thus far, are of interest and are described in the follow- 
ing subsection. 

4.1.1. 3.4 Total nodal mass 

The total mass of a node is given by Equation (10) 



where a = I, r, and c. The masses of inert, reactive, and carbon constituents 
in a node are given by Equations (65), (69), and (70), respectively. As a 
matter of convenience in interpreting computer program output it is appropriate 
that material composition be given in terms of intensive rather than extensive 
properties. The definition of nodal mass fractions and nodal density accom- 
plish this purpose. 




m 


for a = I, r, or c, and the nodal density is given by 


A 6 
n n 


(71) 




4.1.2 Pyrolysis Gas Composition 

In this section, finite difference relations for pyrolysis gas composi- 
tion and flow rate are presented. Most of the fundamental relations required 
for characterizing the pyrolysis gas composition and flow rate were developed 
above in Section 4.1.1, and, as such, the development in this section is rela- 
tively straightforward. 

The gas flow rate leaving a node, n, may be related to the gas flow rate 
entering the node and the gas generation rate in the node by integrating 
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Equation (19) over a nodal volume. 
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From which: 


= m 


am n \ 

W + ^ 


(74) 


It is noted that differentiation of the nodal mass is performed at constant 
y , .and therefore includes only the nodal mass change rate associated with 
chemical reactions. These chemical reactions may consist of gas generation 
resulting from decomposition, or mass transfer resulting from coking reactions. 
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Substituting Equations (47), (48), and (60) into the above yields; 
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where, as before, Ki^ is given by Equations (45) or (46) depending on 
reaction order, and K C g n ) coke is given by either Equations (56) and (57) 
for forward coking,or ^cg n ) co ke = ^cgE n for reve rse coking. 

The elemental carbon content in the pyrolysis gas, K C g* may obtained 
in a number of ways. Perhaps the most straightforward approach is to consider 
the elemental carbon contribution of each of the terms in Equation (76) . The 
gas leaving the node n, has a carbon content, K cg n » and the 9 as entering 
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lias a carbon content K cg n+ q’ The secon< 3 term on the right side of equation 
(76) represents the gas generated in the node as a result of organic decom- 
position reactions. The carbon content of this gas is designated by 

< K c - K c ) 

~ _ i if Lb carbon in gas , ^ 

eg. ~~ \ Lb gas generated 

10 { i j 

c if 


The last term in Equation (76) represents coking. The mass transfer associated 
with coking is pure carbon by definition. Utilizing the above definitions. 
Equation (76) , and performing an elemental carbon balance on the gas results 
in an equation for the elemental carbon content of the gas leaving the node, 
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4.2 CHAR LAYER PRESSURE DROP 

It is desired to cast the pressure drop correlation Equation (29) into 
a form including only variables which are readily available from the computa- 
tion scheme for solving the mass and energy conservation equations. 

The velocity of the gas passing through the porous char is related to 
the gas flow rate per unit area by Equation (30) . Introducing the gas state 
equation (P = p RT/M^) into Equation (30) yields; 
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PM g A 


RT 


(79) 


It is convenient to define a viscosity law for the pyrolysis gas 


U 



(80) 


where 


u q and T q represent a reference viscosity and reference temperature. 
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Substituting Equations (79) and (80) in (29) yields the following rela- 
tion for a linear pressure drop across a node, n. 
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m RT 
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(81) 


where the coefficients a, 3, are functions of the local char permeability. 
The pressure at the top of node 1 (the heated surface node) is taken equal 
to the prescribed boundary layer edge pressure (P^ = P ) . The pressure at 
the top of any other node is related to the pressure drop across and pressure 
at the top of the node above it. 



i=l 


(82) 


(83) 


The computation may be performed in an explicit manner, that is, after obtain- 
ing solutions of the energy and mass conservation equations, a single pass 
through the nodal network from heated surface (n = 1) to the virgin material 
yields the pressure distribution entirely in terms of conditions existing at 
the end of the previous time step. The solution corresponds to assuming an 
impermeable boundary condition at the rear face of the ablative material. The 
summation is performed just to the virgin material since below this point, 
m = 0 . 

g 

4.3 ENERGY CONSERVATION 

A modified differential formulation of the subsurface energy conserva- 
tion equation was developed above, in Section 3.3, and is represented by 
Equation (43) . 
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Finite difference formulation of Equation (43) is relatively straightforward. 
The primary features which distinguish between various possible formulations 
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are associated with whether an implicit or explicit solution is desired. 
Finite difference formulations corresponding to solutions which are explicit 
and implicit with respect to temperature are given here. 


Based upon the satisfactory accuracy checks for the forward difference 
formulation presented in Reference 1, the same basic differencing philosophy 
is adopted here, and results in the following finite difference interpretation 
for Equation (43) . 
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is given by Equation (48) 


is given by Equation (60) 
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is given by Equation (47) 
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The enthalpy of constituent, a, is given by: 



( 86 ) 


The mass fraction is given by Equation (71) 




(71) 


The heat of formation of constituent a, AH f is an input constant and the 

3 , 

specific heat, C , is a function of temperature alone. The pyrolysis gas 

flow rate, m , is given by Equation (76) , and the pyrolysis gas enthalpy, 

hg., is considered to be a function of temperature, pressure, and elemental 

carbon content, K c . Various models may be employed to represent the pyroly- 
g 

sis gas enthalpy, as for example that described earlier in Section 2. What- 
ever model is selected to represent the pyrolysis gas molecular configuration, 
and hence its enthalpy at a given temperature, it would seem that the param- 
eters, p, T, and K c will be required to characterize the pyrolysis gas 

g 

enthalpy. Selection of a model to represent the thermal conductivity in the 
charring material has not been given detailed consideration. Cursory exam- 
ination suggests that the primary factors to be considered in characterizing 
the thermal conductivity would be the relative mass of each constituent, K , 
the porosity, and the temperature. Such a model has been postulated and 
verified to a certain extent (Ref. 9) for partially degraded organic rein- 
forced ablative materials. The same fundamental approach may yield realistic 
correlations for "coked" chars as well. Evaluation of the remainder of the 
terms in the energy Equation (84) is apparent from an examination of the geo- 
metrical configuration being analyzed (Fig. 2) and the nomenclature list. 


Implicit formulation of the energy equation is most conveniently accom- 
plished in the same manner as reported in Reference 1 by linearizing the non- 
linear terms. The following relation is obtained. 
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(Equation continued on following page) 


36 


(Equation continued from preceding page) 
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The specific heats appearing in Equation (87) may be evaluated directly from 
an input tabulation for the material, Cp^, and, for the gas, the specific heat 
is obtained by numerically differentiating the input gas enthalpy data for 
the temperature and pressure at the node, n. Equation (87) is implicit in 
temperature and explicit in mass, that is, the mass transfer events, which 
are temperature dependent, are evaluated at the "old" temperature (at 0) . 

This type of treatment has been found to be quite satisfactory for numerical 
characterization of straight decomposing materials. it is difficult to assess 
the adequacy of this type of formulation for the problem being considered 
here since the kinetic coefficients for characterizing the coking reactions 
are not presently in hand. 

Equations (84) and (87) are appropriate for all nodes in the ablation 
material except the heated surface node, n = 1. Treatment of the heated sur- 
face boundary condition may be accomplished in the same manner as described 
in Reference 1 where it is shown that the surface (char) mass recession rate 
may be related to the surface temperature, surface pressure, pyrolysis off- 
gas rate, and boundary layer mass transfer coefficient. 


B c = f ( P e' T e' B g ) 


where * B g = ™g/p e u e c M A - In the P resent analysis, one addi- 

tional parameter is needed, the pyrolysis gas composition, which may be rep- 
resented by the pyrolysis gas elemental carbon content. 
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Exept for this difference, the surface boundary condition evaluation may be 
treated in the same manner as described in Reference 1. 

SECTION 5 

SUMMARY AND CONCLUSION 

The probable importance of subsurface coking reactions has been identi- 
fied and a phenomenological model to represent energy, mass, and momentum 
transfer events in a charring ablation material undergoing coking reactions 
has been postulated. Differential equations were developed to represent a 
mathematical analog of the phenomenological model and the equations were 
cast into finite difference form suitable for coding into a computational 
scheme. 

The phenomenological model is based upon a certain amount of conjecture 
because sufficient quantitative experimental observations to construct: a. 
rigorous model do not presently exist. Part of the motivation behind conduct- 
ing the study presented herein has been a desire to identify the parameters 
and groups of parameters which require experimental evaluation, since, ideally 
the parameters which are controlled and measured in an experiment may be di- 
rectly related to the parameters required as input to a computational scheme. 
The phenomenological model presented herein is believed reasonable and the 
finite difference representation is believed appropriate for developing an 
efficient computational scheme. 
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Initial pyrolysis gas composition 
(C,H,0, „ = 0.62039) 



Fiqure 1. Gas Phase Elemental Carbon Concentration and Coking 
Potential for Typical Phenolic Pyrolysis Gas 
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Figure 3. Experimentally Determined Coefficients for 
Flow through Porous Media 
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APPENDIX A 


DEVELOPMENT OF RATE EQUATION FOR CHARACTERIZING 
DECOMPOSITION OF ORGANIC CONSTITUENTS 

The decomposition rates of organic constituents utilized in ablative 
materials have been successfully correlated in terms of kinetic equations of 
the Arrhenius form (Refs. A-l and A- 2) „ The purpose of this appendix is to 
develop a kinetic equation of a form convenient for inclusion in a charring 
material response calculation procedure in terms of quantities readily de- 
rived from TGA data. 

TGA data are normally interpreted in terms of an idealized irreversible 
reaction of the form 

k f 

m — ► m + m (A-l) 

eg v 1 

where m represents the instantaneous sample weight which is decomposing to 
form a carbonaceous residue, m c , and pyrolysis gas, m g . The TGA data are 
correlated by an equation of the form 

d / m , -E/RT f m 

air hr = - k f e hr 

l o ^ o 

where 

m = instantaneous sample weight 

m Q = initial sample weight 

m cf = weight of carbonaceous residue after decomposit ion 

= pre-exponential constant for the reaction 
n = reaction order 

E = activation energy 

The reaction rate coefficients in the above equation are derived directly 
from TGA data, however, the equation is not suitable, in its present form, 
for inclusion in a computational scheme where consideration is given to a 
number of simultaneous decomposition reactions. This is the case because 
Eq. (A- 2) contains the carbon residue mass (m . If decomposition of a 

mixture of constituents is considered, each characterized by Eq. (A- 2) , thou 
it would be necessary to keep track of the carbon residue that has resulted 
from decomposition of each constituent. This requires extra computational 
effort, which, as will be shown, is not necessary. It is desired to express 
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the decomposition rate of each reactant in terms of the reactant mass rather 
than the total mass of reactant and char residue. In Eq. (A-2) the total 
mass at any instant, m, represents the mass of undecomposed reactant, m^, 
plus the mass of char residue, m c * 

m = m + m (A- 3) 

r c 

The sketch below depicts the masses m and m as a function of the total 

r c 

mass, m 



It is noted that initially, before decomposition. 


m = m r = m Q 

and finally, after decomposition is complete 


m = rn ^ = in 

cf c 

The rate of consumption of reactant may be related to the total mass loss 
rate as follows 
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Utilizing the above sketch it is 
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Substituting (A-4) in (A-5) yields 
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Substituting (A-6) in (A- 2) results in an expression for the rate of reactant 
consumption 
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Rearranging and integrating Eq. ( A— 5) obtains 
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Substituting in Eq. (A-7) yields 
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The above is the desired form of the decomposition rate equation for a system 
having several reactive constituents because it contains neither the total 
mass, m, nor the mass of char residue at any instant, m . Utilization of 
Eq. (A- 2) directly would require keeping track of the quantity of carbon 
residue resulting from each of a number of decomposition reactions. Since 
the decomposition rate clearly depends upon the quantity of reactant present 
and not upon the quantity of carbon previously formed from the reaction, 

Eq. (A- 10) is simpler to use in a computational scheme. 

Rewriting Eq. (A-10) for the i th reactant and defining mass fractions 
on the basis of the initial mass, m Q , yields 
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(A— 11) 


where 
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mass of i per unit mass of i prior to 
decomposition 

(m /m. ) mass of carbon residue after complete 
c^ 10 

decomposition per unit mass of i prior to decomposition 


It is noted that both k . and K are constant for a given reactant, 

ri c if 

and # as such, it is convenient to combine them into a single effective rate 
coefficient, k^. 
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